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Abstract 

The formation of a thin current sheet in a magnetic quasi-separatrix layer (QSL) is investi- 
gated by means of numerical simulation using a simplified ideal, low-/?, MHD model. The initial 
configuration and driving boundary conditions are relevant to phenomena observed in the solar 
corona and were studied earlier by Aulanier et al., A&A 444, 961 (2005). In extension to 
that work, we use the technique of adaptive mesh refinement (AMR) to significantly enhance the 
local spatial resolution of the current sheet during its formation, which enables us to follow the 
evolution into a later stage. Our simulations are in good agreement with the results of Aulanier 
et al. up to the calculated time in that work. In a later phase, we observe a basically unarrested 
collapse of the sheet to length scales that are more than one order of magnitude smaller than those 
reported earlier. The current density attains correspondingly larger maximum values within the 
sheet. During this thinning process, which is finally limited by lack of resolution even in the AMR 
studies, the current sheet moves upward, following a global expansion of the magnetic structure 
during the quasi-static evolution. The sheet is locally one-dimensional and the plasma flow in its 
vicinity when transformed into a co-moving frame, qualitatively resembles a stagnation point flow. 
In conclusion, our simulations support the idea that extremely high current densities are generated 
in the vicinities of QSLs as a response to external perturbations, with no sign of saturation. 
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I. INTRODUCTION 



The spontaneous formation of thin current sheets is believed to play an important role 
in the dynamics of astrophysical plasmas like the solar corona and, in particular, for the 
onset of magnetic reconnection. In fact, impulsive events like solar flares and coronal mass 
ejections are often associated with magnetic reconnection as a mechanism that effectively 
releases magnetic energy, which in turn calls for highly concentrated electrical currents to 
explain the breakdown of ideal plasma behavior by means of e.g. micro-instabilities in the 
non-collisional coronal plasma. 

Recently, the study of magnetic field configurations as candidates for reconnection has 
shifted from separatrix surfaces towards quasi-separatrix layers (QSLs) [l-3]. In contrast 
to genuine separatrix configurations, QSL fields do not necessarily contain magnetic null 
points in 3D, which makes them relevant in many more situations than strict separatrices. 
QSLs describe a magnetic field mapping between two boundaries that is continuous, but 
changes rapidly in space. This change has been quantified in terms of a flux tube "squashing 
factor" Q 4] and more recently generalized to remove boundary projection effects [5|. In 
accordance with the significance of QSLs for magnetic reconnection, the problem of current 

n n 

sheet formation has been investigated. Here, Titov et al. [6j and Galsgaard et al. [7J studied 
the formation of current sheets in a straight hyperbolic flux tube (HFT) configuration both 
analytically and numerically, and found exponential growth of the current density as a 
response to shearing magnetic footpoint motion. One major conclusion of that work was 
that a shear in the applied boundary motion was important, if not essential, for the formation 
of thin current sheets, and that the creation of a stagnation point flow in the HFT was the 

cey element in this process. However, this has been later questioned by Aulanier et al. 

8): They argued that the initial squashing factor Q in the simulations of Galsgaard et al., 
probably together with additional symmetry properties, was too small to account for highly 
concentrated currents, and that the shear boundary motion and stagnation point flow in 
[3] served as to dynamically create thin QSLs only during the simulation. This argument 
was underpinned in jg] through extensive comparison between MHD simulations of less 
symmetric potential magnetic field configurations which initially result from four magnetic 
point sources submerged below a photospheric boundary. They contain QSLs with squashing 
factors of up to Q ~ 10 5 and were exposed to different boundary driver patterns, namely 
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a shear and a translational motion. Aulanier et al. observed the formation of thin, intense 
current sheets, almost irrespective of the type of applied boundary driving, which stresses 
the relevance of the initial field geometry and QSL strength, rather than the boundary 
motion, for current sheet formation. Later, that work has been extended by finite resistivity 



to simulate magnetic reconnection at that thin current sheet jsj with its strong temporal 
change in magnetic connectivity. 

A general limitation of these previous numerical studies was the lack of reasonable numer- 



ical grid resolution at the sites of current sheet formation [9] : During the formation process, 
the currents get more and more concentrated locally so that the sheet scales quickly reach 
the numerical grid spacing. This left open the question whether the current concentration 
would continue to length scales small enough to account for the onset of micro-instabilities, 
or if this process would get arrested at some larger length scale. In order to address this 
question, but also get more insight into the local dynamics in the immediate vicinities of 
the current sheets, we have carried out numerical simulations in the spirit of js], using one 
of their initial configurations and boundary drivers. However, we used the technique of 
adaptive mesh refinement (AMR), implemented in the code racoon jlp| . in order to obtain 
a much higher grid resolution in the current sheet vicinity than was reported earlier. This 
allowed us to obtain insight into the formation process that reaches beyond the previous 
results. 

In the next section, we will briefly describe the model that we employed in our studies, 



while results that complement the work of [8] are given in Sec. Ill, followed by a summary 
and conclusion in the last section. 

II. BASIC EQUATIONS AND NUMERICAL MODEL 

Our simulations are based on a reduced subset of the MHD equations, appropriate for 
the quasi-static, low-/3 regime, where /3 is the ratio of thermal to magnetic pressure. In 
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normalized quantities, it reads 

d t v = —v ■ Vv + —j x B + uAv (la) 

d t B = V x (v x B) + r?AB - Vtf (lb) 

j = V x B (lc) 

= -c^V-B + qA* (Id) 

Pressure terms are omitted from the momentum equation ffTal) because of the low-/3 ap- 
proximation. As we are interested in a quasi-stationary evolution of the system, i.e. the 
limit of high wave speeds, we replace the continuity equation by a direct prescription of a 
relaxation mass density, namely p := B 2 . This approach results in an homogeneous Alfven 
velocity ca '■= \B\/y/p = 1 and fast communication of unbalanced forces throughout the 



system, which shortens relaxation times and has been successfully used in other studies [11]. 
Constant kinematic viscosity v — 5 • 10~ 4 and resistivity rj = 5 • 10~ 6 are included only to 
guarantee numerical stability on the grid scale and have little effect on the overall evolution. 
The artificial scalar function \l/ and its dynamic equation ( lldj) serve as a convenient means 



to constrain any finite V • B, resulting from discretization errors, to negligible values 
Combining <9 t V • fllbj) and A (lldj) results in the mixed equation 



12|: 



d tt V ■ B = c|AV • B + (r) + ci)Ad t V ■ B - q//A 2 V ■ B 

for V • B. The crucial term on the right side is the first Laplacian, which describes a 
hyperbolic transport of V • B with velocity Ch and leads to radiative distribution of V ■ B 
throughout the computational domain, while it gets dissipated by phase mixing and diffusion 
according to the other two terms. Note that there is freedom in a specific choice for the 
\l/-equation flld|) a s it is of the order of the discretization error anyway. For instance, while 
Dedner et al. [12] discuss the term — %^ (see their Eq. (19) resp. (24e)) to arrive at 
a telegraph equation for \l/ and V • B in the continuous case, we found this unnecessary, 
although possible, for our computations. The main reason for this seems to be that in our 
case, the sources of V • S-errors are highly localized regions in space, namely the regions of 
intense currents, so that the hyperbolic transport is the dominating cleaning effect. On the 
other hand, we added the term qA\& in Eq. (Ildl) . which is not discussed in this specific form 
in 12|. Its motivation is, however, not so much a change in the V • -B-cleaning property 



itself, but the observation that the purely hyperbolic choice, i.e. Eq. fjldl) with q = 0, tends 



to introduce odd-even-decoupling on the centered finite difference grid that we used. This 



decoupling, which was not an issue in 12| since they employed finite volume schemes in finite 
element discretisations, could be healed satisfactorily through the additional Laplacian term 
which couples values on neighboring grid points with each other. We finally found overall 



good V • .B-cleaning properties when choosing the numerical parameter values to q = 5 • 10 
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and 4 = 5- 1CT 2 . 

The equations are discretized in a domain of (x,y,z) G [—0.7,0.7] x [—0.5,0.5] x [0,0.5], 
where we identify the plane z = with the solar photosphere. Integration is performed with 
a strongly stable third-order Runge-Kutta scheme [13j , spatial differentiation is realized with 
standard second-order finite differences. In order to resolve the expected small-scale features, 
we employed the block-structured adaptive mesh refinement (AMR) framework racoon jlo| . 
Here, we used the norm of the magnetic field gradient, VB, to derive a local length scale 
for B that serves as an indicator criterion for local mesh refinement. Effective local grid 
resolutions obtained with this technique were 4096 3 in the present work. 

Initial conditions are adapted from [8], where we address the magnetic field configuration 
resulting from two pairs of opposite polarity photospheric flux patches whose respective 
connecting axes intersect at an angle of 150 degrees. Specifically, the initial magnetic field 
stems from 4 virtual magnetic point sources, indexed by i, below the photosphere 

i=l 

with respective source strengths F\ = —F 2 = 1 and F 3 = — F 4 = 0.4 and locations 
7*1,2 = 0, — |) and r 3 4 = (T^, — th)j respectively. This field geometry is known to 
contain two symmetric dome-shaped QSLs with squashing factors of Q ~ 10 5 , intersecting 
in a hyperbolic flux tube (see js| for details). 

In the course of the simulation, the field is exposed to a horizontal photospheric vortex 
flow around the magnetic source % = 3 in Eq. (j2J), realized by prescribing the boundary 
condition for v at z = 0. The maximum flow velocity is msix(\vBnd.\) ~ 2 • 10 -2 and it gets 
ramped up in time according to oc ~ [tanh(|(t — 1)) + l]. 

The detailed treatment of the lower boundary is as follows: As all quantities are dis- 
cretized as cell-centered variables, boundary conditions have to provide values for v, B and 
^ at Z\/2 '■= A 2 /2 in a way that is consistent with the evolution equations and the overall 
second order accuracy in the grid spacing A z . Denoting the boundary values of individual 



quantities at z = with a hat, the velocity v = VBnd. is explicitly given from the prescribed 
horizontal vortex flow, and in particular, v z = 0. All components of v at zui can now 
simply be interpolated in ^-direction with second order accuracy between z = and the 
interior grid points above Zi/2- Using this direct forcing, there is no need to compute j or 
the Lorentz force at z\ii at all, hence we do not need the horizontal components of B at this 
stage. Now, to integrate B, we note that the convective part of the _B 2 -equation becomes 
autonomous in the boundary plane, involving no other undetermined quantities nor any 
derivatives normal to the boundary: Writing the convective electric field as £ := — v x B 
gives S x = —v y B z and £ y = v x B z so that d t B z = — ■ (vhB z ) where the index h indicates 
horizontal vector projections, e.g. the prescribed Vh = (v x ,v y ), and := (d x ,d y ). In other 
words, when ignoring the numerically motivated resistive diffusion and V^-terms, we are 
able to integrate the "proper" B z completely for its own, which in turn allows us again to 
interpolate B z at Z\/2 up to 0(A 2 Z ), similar to v above. Now, this "proper" £? z -equation also 
determines £h, which allows to update Bh at z\/i- At this point, the magnetic diffusion 
and V ■ -B-cleaning terms are taken into account as usual, where the former requires an 
extrapolation of Bh across z = 0. Finally, the right side of Eq. ( Ildl) is easily evaluated at 
Zi/2, using B z from above and the Dirichlet condition \1/ = 0. 

While the upper and lateral boundaries can in principle be handled in the same fashion, 
using the no-slip and no-penetration condition v = 0, we used a simpler approach there: 
Setting v to zero ahead of the boundaries, keeping the tangential components of B fixed 
in ghost cells and applying solenoidal and homogeneous Dirichlet conditions to the normal 
components of B and respectively, proved to be sufficient for those passive boundaries. 
Note, in particular, that no artificial Lorentz forces act there either. 

The boundary treatment described above presumes a homogeneous numerical grid. It 
has been applied in the same spirit to the AMR simulations that we present here, where 
additional complications occur at the interfaces between neighboring grid blocks of different 
mesh resolution that abut the physical boundaries. Without going into details, we only 
mention that additional coarse-fine and fine-coarse interpolations are needed for the com- 
putation of B z and £ at those junctions, but that they do not pose any fundamental new 
challenge apart from the programming complexity. 
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FIG. 1: Color coded \j\ in the plane y = at t = 5.0. Extended current systems match those 
reported by 8|] . The thin current sheet of interest formed inside of the red rectangle and is hardly 
visible on these scales. Maximum values of \j\ are ~ 400 near the photospheric boundary, and 
~ 1000 in the marked region. Block lines indicate the layout of the recursively refined grid blocks, 
each containing 16 3 cells. The two finest block levels are omitted for clarity. 



III. RESULTS 



After applying the photospheric boundary driving in v, dynamic shear modes travel into 
the domain and mix there. On the scale of a few Alfven transit times, the magnitude of the 
current density grows significantly and a quasi-stationary current system as shown in Fig. [T] 
builds up. It consists partly of relatively weak currents which are distributed on a large 
scale in a dome-like structure that is pre-determined by magnetic field lines connecting the 
driver region with the opposite polarity regions of the photosphere. On top of this, a highly 
localized thin current sheet can be identified in the vicinity of (x,z) w (0,0.18) in Fig. [IJ 
This thin current sheet actually lies inside the pre-existing QSL of the initial magnetic field. 
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FIG. 2: Plot of \j\ vs. z at x = y = taken at times t = 4.2,4.7,5.0,5.4,5.8. The inset expands 
the relevant current sheet height range. Maximum values of |j| are ~ 300, 600, 960, 1700 and 2900, 
respectively, increasing monotonically in time. The current sheet moves upwards and thins with 
respective FWHM values of « 4.2 • 10~ 3 , 1.9 • 10~ 3 , 1.2 • 10" 3 , 7.2 • 10~ 4 and 4.8 • 10" 4 . 



We interpret the striation patterns at x ~ —0.4 in Fig. [T] as signatures of MHD waves on 
field lines which connect the strong photospheric field region with the current sheet during 
the evolution. 

These results are basically in good agreement with those published earlier by Aulanier et 
al. 8j. However, the sheet thickness at the stage shown in Fig. [1] is already on the numerical 
grid scale of js], so that their studies were unable to investigate into the further evolution 
or features on smaller scales. 

The temporal evolution of the sheet is displayed in Fig. |2] which shows the vertical profiles 
of \j\ at x = y = for four different times. It is evident that the sheet thickness decreases 
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FIG. 3: Growth of the maximum of \j\, taken over the entire domain, against time. 

with time while the value of maximum current density increases accordingly to reach values 
of ~ 3 ■ 10 3 in the latest stage. At the same time, the sheet moves upwards, as indicated 
by the inset graphs. In fact, we find that the entire magnetic structure expands gradually 
as a response to the boundary driving, and the current sheet as a substructure is embedded 
in this motion. An other detail that emerges from Fig. [2] is the fact that up to t w 4.4, the 
most intense currents are not yet found in the thin current sheet itself, but in the large-scale 
system at z ~ 0.02, i.e. close to the photospheric driver (see also Fig. [1]). It is only at later 
times, that the thin sheet dominates in the current intensity. 

This phenomenon is also visible in the temporal evolution of max(|j|), which is plotted 
logarithmically against time in Fig. |3j The early phase, with rather fast growth of maxdj'l), 
corresponds to the ramp-up of the boundary driver, which essentially reaches its maximum 
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FIG. 4: Color coded \j\ and velocity components (v x ,v z ) as arrows in the plane y = at times 
t = 4.2,5.0 and 5.8 (left to right). Arrows in the upper row show the plasma velocity in the 
fixed reference frame, while v has been transformed into the co-moving frames of the current sheet 
for the lower figures. The transformation velocities are (V X ,V Z ) = (—3,6) • 10 _3 ,(— 10,8) • 10 -3 
and (—7,7) • 10 -3 , respectively with the transformed |(ur,i>z)| attaining maximum values of 1.6 • 
10~ 2 ,2.3 • 10~ 2 and 3.0 • 10~ 2 (left to right). Note also that the x- and z-coordinates on the axes 
are relative to the point (xq, zq) = (—2.5, 18) • 10~ 2 . 



magnitude around t w 1.7. This is followed by a slower growth rate of the current maximum 
up to t m 4.7. During this stage, the maximum values stem from the extended currents 
close to the photosphere (compare with Fig. [2]), which eventually get overtaken by the 
faster growth of the thin embedded current sheet. Further intensification continues, with 
amplification of max(|j|) by roughly a factor of 5, until the growth slows down at t ~ 5.5. 
At this time, the sheet thickness is only a few times the numerical grid scale and thereby 
poorly resolved with artificial diffusion effects becoming competitive. 

Fig. H] shows details of the sheet in the cut plane y = for three different times. Again, 
the overall upward motion and the thinning and intensification of the sheet are well visible. 
In addition, we have plotted the plasma velocity as arrows, projected into the displayed 
plane, to give an impression of the flow in the vicinity of the current sheet. While the upper 



10 



row shows the velocity relative to the fixed simulation frame, the plots in the lower row 
show the flow transformed to a frame which is co-moving with the expansion velocity of the 
structure. To this end, we identified the locations of maximum current density in the plane 
y = from plots of successive data output sets, and computed a pattern velocity from their 
difference. This velocity was then subtracted from the plasma flow velocity in the lower 
plots of Fig. 4. There have been controversial discussions as to whether the current sheet 
formation at the hyperbolic flux tube embedded in the QSL is related to a stagnation-point 
flow 7|, |8j. In particular, Aulanier et al. claim that no stagnation point exists at the intense 
current sheet. This is certainly confirmed by our simulation, however we remark that the 
focus on a strict definition of stagnation point, i.e. v = 0, maybe somewhat misleading 
because i) the velocity is sensitive to the chosen frame of reference, and ii) the component 
along the current sheet should be discarded from these considerations anyway, because it 
largely decouples from the mechanism of magnetic shearing discussed in [6| and The 
flow pattern projected into the y-z-pl&ne is shown in Fig. [5j again transformed into a frame 
that moves upward with the current sheet and, in addition, results in v y = in the current 
sheet center. This figure also demonstrates that the current sheet is indeed elongated in the 
y-direction. Hence, at least in the latest stage displayed in Figs. H] and O it can be treated as 
a quasi-one- dimensional sheet. Finally, we remark that the assumption of a quasi- stationary 
evolution loses its validity in the late stage of the sheet evolution: Obviously, the collapse 
becomes a localized, dynamic process associated with significant magnetic forces. This can 
be seen from the field line plot shown in Fig. [6J where magnetic field lines connecting the 
thin current sheet with the photosphere have been colored with the quantity a := 
For a force- free field, j x B = 0, the value of a is constant along magnetic field lines. This 
condition is obviously not met in the QSL, which means that the currents close locally across 
field lines. 



IV. CONCLUSIONS 



We carried out numerical simulations of current sheet formation in a quasi-separatrix layer 

using a simplified MHD model appropriate for the quasi-static evolution of a low-/? plasma. 

I I 

The setting under consideration has been investigated before [8] and our results agree well 
with that work as long as the current sheet structure is well resolved in both studies. With 
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FIG. 5: Color coded \j\ and velocity components (v y ,v z ) at x = —2.5 • 10 -2 and t = 5.8, cor- 
responding to the bottom right plot in Fig. HJ Here, v has been transformed into the co-moving 
frame (V y , V z ) = (—6, 7) • 10 -3 and max(|v|) ~ 1.1 • 10 -2 in that frame. Coordinates are relative to 
{y ,z ) = (0,0.18). 

the use of local adaptive mesh refinement (AMR), we were able to follow the thinning of 
the current sheet further down to a scale which is about one order of magnitude smaller 
than previously investigated. In particular, our simulations reached a stage in which the 
maxima of \j\ in the upper part of the QSL grow significantly beyond the values close to the 
photospheric boundaries, which gives clear evidence that the most intense current densities 
actually can be expected in the upper part of the QSL. This late stage is characterized by 
a relatively fast collapse of the locally almost one-dimensional sheet with an approximately 
exponential increase of max(|j|) in time, and the evolution is no more quasi-static at this 
point. 

When magnetic forces become significant in this late stage of the current sheet formation, 
full nonlinear MHD dynamics will take place. Previous studies have addressed details of the 
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FIG. 6: Color coded a = \j\/\B\ in the planes y = and z = at t = 5.8. The maximum 
value a m ~ 2 • 10 3 is attained in the current sheet center (red). The magnetic field lines, starting 
equidistant from (—0.03,0,0.175) to (—0.03,0,0.185), are also color coded with a and show that 
B ■ Va 7^ 0, i.e. the magnetic field deviates significantly from a force-free field. 



local dynamics of such current sheets using appropriate initial conditions and periodic sys- 
tems (e.g. [ijl and references therein). There, one particular question has been whether the 
current density might form a singularity in finite time, or whether its growth is limited to 
merely e.g. exponential behavior. On theoretical grounds, it could be shown that a dynam- 
ical alignment between the velocity and the magnetic field would bound \j\ to exponential 
growth. Analyzing our data further, we actually found indications of such an alignment (not 
shown here), so that we expect to see a collapse of the sheet with exponential growth, i.e. a 
continuation of the phase observed between t m 4.7 and 5.5 in Fig. [31 given that it could 
be resolved numerically. At present however, even our AMR simulations are limited by the 
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lack of further resolution and by numerical side-effects like artificial stabilizing diffusion. 

The plasma flow pattern has been analyzed in a cut plane that is approximately perpen- 
dicular to the local direction of the current density in the sheet: If transformed into a frame 
that moves with the expanding structure, it exhibits a shear pattern which arises from a 
vortex below and a large-scale flow above. The vortex maps to the vortical boundary driver, 
while the large-scale flow is related to the slow overall expansion of the magnetic structure. 

In this paper, we have only addressed the case of one specific boundary perturbation, 
namely a twisting motion at the footpoints. We have also conducted a number of simulations 
with a translational motion analogous to that used in |8j , and observed comparable behavior 
in these well. In particular, the achieved current densities continued to rise at similar 

rates until they were restricted by finite grid resolution. 
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